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Proposed space-based gravitational wave antennas involve satellites arrayed either in an equilat- 
eral triangle around the earth in the ecliptic plane (the ecliptic-plane option) or in an equilateral 
triangle orbiting the sun in such a way that the plane of the triangle is tilted at 60° relative to the 
ecliptic (the precessing-plane option). In this paper, we explore the angular resolution of these two 
classes of detectors for two kinds of sources (essentially monochromatic compact binaries and co- 
alescing massive-black-hole binaries) using time-domain expressions for the gravitational waveform 
'(«--N ' that are accurate to 4/2 PN order. Our results display an interesting effect not previously reported 

f^ , in the literature, and particularly underline the importance of including the higher-order PN terms 

f^ ■ in the waveform when predicting the angular resolution of ecliptic-plane detector arrays. 
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I. INTRODUCTION 

The era of gravitational astronomy is nearly here. 
In the next few years, ground-based laser gravitational 
wave detectors will begin to come on line and, it is 
hoped, space-based detectors will soon follow. Even if 
the ground-based detectors do not succeed in detecting 
the waves within the next decade, the space-based detec- 
tors will certainly be able to, so the physics goal of the 
detection of the new phenomenon will be achieved one 
way or the other, and gravitational wave detectors will 
then move from being physics experiments to being tools 
of observational astronomy. 

It can be argued that the frequency band accessible 
to the space detectors is the more interesting one, both 
experimentally and observationally. It is here, at frequen- 
cies between 10^^ Hz and 1 Hz, that identified individual 
sources (such as AM CVn) exist, thus guaranteeing the 
detection of the phenomenon. But it is also at these fre- 
quencies that astronomical objects of the greatest obser- 
vational interest are to be found, objects for which only 
gravitational astronomy can provide data that bears on 
critical astrophysical and cosmological questions. In this 
paper, we will investigate the ability of space-based de- 
tectors to determine the locations in the sky of two of 
these kinds of astronomical sources — close compact bi- 
naries in the Galaxy and coalescing massive black hole 
binaries in the nuclei of galaxies at high redshift. 



A. Gravitational Astronomy Science Goals 

There are important questions relative to the statistics 
and distribution of close compact binaries in the Galaxy. 
Because of their low luminosity, only a few nearby white 
dwarf binaries have been observed, and even fewer neu- 
tron star binaries have been observed, due mainly to their 
narrow beams as radio pulsars. However, these endpoints 
of massive close binary evolution are the touchstone for 
the binary evolutionary models. Good average density 
statistics are needed in order to determine parameters 
of the models and to provide the basis for calculation of 
the branching ratios between close binary formation and 
other possible events such as type I supernovae. Space 
gravitational wave detectors are omnidirectional instru- 
ments that will provide a complete survey of close bi- 
naries. However, determining the statistics depends on 
measuring the spatial density, which is determined by 
solving for the positions of the objects that are seen grav- 
itationally. The ability of a space detector to determine 
the position of the binary system it is observing is there- 
fore an important aspect of the scientific value of the 
detector. 

Although more speculative than the close compact bi- 
naries, there is another source that is expected to be even 
more interesting for physics, astronomy, and cosmology. 
Evidence continues to mount for the existence of massive 



black holes in the nuclei of most galaxies, including our 
own, and there is increasing evidence that most galaxies 
have been involved at some time in their past in a col- 
lision and merger with another galaxy. If both of these 
things are true, then an exciting source of gravitational 
waves is predicted — the coalescence of a black hole bi- 
nary, formed by the collision of galaxies and the subse- 
quent sinking of the black holes to a common center via 
dynamical friction. The event rate to be expected is a 
few per year for events out to z = 3. These sources will 
be so bright, even at this distance, that they will domi- 
nate the noise spectrum in the sensitivity window of the 
instruments and will be seen for months prior to the fi- 
nal coalescence event. In addition, if the massive black 
holes themselves initially grew by many small mergers 
of this type, then event rates of several tens of visible 
events per year would be expected. Observations of coa- 
lescences will provide valuable physics, since the observed 
waveforms represent ground truth for the computer codes 
presently being constructed to model the dynamics of 
these strong-field events. They will also provide insight 
into galactic merger rates, galactic and protogalactic cos- 
mogeny, etc.. 

Much of the scientific value from the detection of such 
an event depends simply on observing the signal and 
determining from the waveform what the parameters of 
the source are, regardless of where the source is located. 
However, an important aspect of the science of the coa- 
lescences could also be to detect an optical counterpart 
of the galaxy in which the event is taking place, to see 
if there are observable effects produced in the galaxy by 
such an energetic event in its nucleus. For this goal, a 
position for the source would need to be determined. As 
we shall see, the predicted angular resolution of the space 
detectors will probably not be good enough to identify a 
single galaxy associated with the coalescence. However, 
it should be possible to identify the correct portion of the 
sky well before the final event, so that the region may be 
monitored and the galaxy identified by looking for some 
other anomoly in it. Finally, if there are many events 
seen per year, a knowledge of the spatial distribution of 
coalescences would be of value in understanding how they 
relate to galaxy formation statistics. 

The purpose of this paper is to investigate the posi- 
tional sensitivity of space gravitational wave detectors 
when observing signals from these two classes of sources. 
With signals at the expected levels for realistic sources 
and with noise at the levels expected for the instruments, 
we perform a least-squares covariance study for the pa- 
rameters representing the direction to the source. 



B. Space-Based Gravitational Wave Detectors 

The geometry of the detectors is shown in Figure W. 
Three spacecraft|^ foUow roughly circular orbits at the 
vertices of an equilateral triangle inscribed on the circle. 
The spacecraft at each of the vertices can function as 
the central mirror of a Michelson interferometer, while 
the two far spacecraft that it tracks function as the end 
mirrors of such an instrument. The passage of a grav- 
itational wave changes the curvature of space between 
the three spacecraft and affects the phase of the laser 
signals passing between them, producing a signal in the 
interferometer. Two of the three interferometers are lin- 
early independent of each other and provide independent 
information on the gravitational waveform. 

The response of these detectors to a passing gravita- 
tional wave depends on the orientation of the interfer- 
ometer relative to the source. We consider two different 
orientation scenarios. In the first, the plane of the or- 
bit is fixed in space and the triangular geometry simply 
rotates in the plane at some orbital period. This is the 
geometry of the proposed OMEGA mission in which the 
spacecraft orbit the earth in a plane close to the ecliptic. 
In the second scenario, the plane of the orbit precesses 
and the triangle rotates in the precessing plane with the 
same period. This is the geometry that has been studied 
for the LISA mission, acheived by placing each space- 
craft in a heliocentric orbit that is slightly eccentric and 
slightly inclined, so that the spacecraft all remain in a 
plane that is inclined 60° to the ecliptic. The dynamics 
of the heliocentric orbits then cause this plane to precess 
about the normal to the ecliptic once per year. Both or- 
bits are sensitive to the position of the source as a result 
of their motion. 



II. DIRECTIONAL SENSITIVITY 

A. Extracting Directional Information 

It is very difficult to determine the direction to a source 
that emits only a short burst of gravitational waves. 
Since the magnitude of the detected signal depends on 
both the amplitude of the wave and its direction of prop- 
agation relative to the detector arms, one would need to 
know the intrinsic amplitude of the wave (which would 
involve, for example, knowing accurately the distance 



to the source) to say much about its direction. Even 
then, the direction of the source is underdetermined, as 
there are many different propagation directions that yield 
the same signal strength. However, the orbiting-binary 
sources considered in this paper produce gravitational 
wave trains that have predictable (slowly changing) fre- 
quency and amplitude and undergo many oscillations 
during the course of a year of observation. Over this time 
span, two distinct effects modulate the detected signal in 
a manner that depends on a source's position in the sky. 

The first of these effects is the Doppler shift due to the 
detector's orbital motion around the sun. The phase of 
this Doppler shift depends on the azimuth angle of the 
source in the plane of the ecliptic (the azimuthal direction 
of the source, that is, the direction of the projection on 
the ecliptic plane of a vector pointing toward the source, 
is the direction that the detector is moving when the 
source is maximally blue-shifted) . The amplitude of this 
yearly Doppler shift depends on the polar angle of the 
source relative to the zenith of the ecliptic (clearly, if the 
source were at the zenith, there would be no Doppler shift 
at all, while if the source were in the plane of the ecliptic, 
the amplitude of the Doppler shift would be maximal). 
So measuring the phase and amplitude of the observed 
Doppler shift of the wave train allows one to determine 
(to some accuracy) the angular coordinates of the source. 

If the detectors are deployed in the precessing-plane 
configuration, then a second effect comes into play. As 
previously discussed, the normal vector to the plane of 
the detector array in the precessing-plane case sweeps 
around the sky as the array orbits the sun. Since the de- 
tected signal depends on the direction of the wave prop- 
agation relative to the detector array, its amplitude will 
thus be modulated in a distinctive way as the plane pre- 
cesses. Monitoring the change in signal amplitude as 
a function of time thus gives us additional information 
about the angular coordinates of the source. Since the 
normal to the detector plane is fixed in the ecliptic-plane 
case, one gets no information from this effect in that case. 

The quality of the information that can be extracted 
from the Doppler-shift effect improves as the frequency 
of the source increases (having a larger number of wave 
cycles during a given time improves one's ability to deter- 
mine the Doppler shift). The quality of the information 
obtained from the second effect, by contrast, is compar- 
atively insensitive to frequency (as long as the wave still 
undergoes many cycles per year). We will see in Sec- 



tion 



II] that our results are consistent with these general 



principles. 



^The currently proposed LISA mission presently has just 
such a configuration, in which each of the three spacecraft 
must simultaneously track in two directions. The OMEGA 
mission and a previous configuration of the LISA mission use 
two separate spacecraft at each vertex of the geometry, so 
that each spacecraft can rigidly orient itself to track a single 
direction. 



B. Goals for This Project 

The question of how well one can determine the an- 
gular coordinates of binary sources using these effects 
has been previously examined by Cutler ||] and Cutler 
and Vecchio PS]. Cutler's first paper discusses the an- 



gular resolution of precessing-plane (LISA) option only, 
while the papers by Cutler and Vecchio mention briefly 
the ecliptic-plane (OMEGA) option as well. In all of 
these papers, the gravitational wave is analyzed in the 
frequency regime, and while the wave phase is modeled 
to 3/2 post-Newtonian order (in the first paper) or to 
4/2 PN order (in the latter two), the gravitational wave- 
form itself is modeled only to lowest PN order (that is, 
it includes only a single harmonic). Also in these pre- 
vious papers, the parameter space associated with the 
angular parameters is sampled randomly, and (in dealing 
with very distant events) the universe was assumed to be 
empty for the sake of simplicity. 

In the work described here, we sought to extend and 
improve on this previous work by 

1. Extending the approximation to 4/2-PN order for 
both the phase and the waveform, 



2. Treating the universe as being flat {fl 
instead of empty (fi = A = 0), 



l,A-0) 



3. Systematically exploring the parameter space to ex- 
pose patterns, and 

4. Comparing the precessing-plane and the ecliptic- 
plane cases more fully. 

We also opted to work in the time domain instead of 
the frequency regime and calculate all derivatives ana- 
lytically instead of numerically. 

Our effort to include higher-order terms in the wave- 
form was driven by a desire to more accurately model the 
coalescing black-hole scenario, and it turns out (as dis- 
cussed in section M below) that including these terms is 
especially important in the ecliptic-plane case. Although 
it is more complicated, we considered the flat, Friedmann 
universe to be a better model of the real universe than 
the empty universe (the distinction might be important 
for handling black-hole coalescence events at very large 
redshifts). We also felt that systematically exploring the 
parameter space would help us understand the behavior 
of these two detector conflgurations more physically. 



III. METHOD 

A. The Received Gravitational Wave to 4/2 Order 

Blanchet, et. al. B provide a complete description of 
the gravitational waves emitted by inspiraling compact 
binaries to 4/2-PN order, assuming that the objects have 
zero spin, are essentially point masses, and are in a qua- 
sicircular orbit. The waveform as observed in the frame 
of the detector array can be expressed as follows: 
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where r is the geometrical distance to the source, c is the 
speed of hght. 
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expresses the system's total mass in units of time, 
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is a unitless expression of the ratio of the objects' masses 
(?7 = 1/4 when the masses are equal but 77 — > as one of 
the masses becomes very large compared to the other), 
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(where oj is the time-dependent angular frequency of the 
binary's orbit in its own frame) is a time-dependent ex- 
pansion parameter that is of the same order of magnitude 
as the ratio v /c for the system, and 
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The complete expressions for the higher-o rder Hs are 
found in Blanchet, et. al. [Q. In equations 3.3, i is the 
angle of the source's orbital angular momentum vector 
relative to a unit vector pointing from the source to the 
detector (z is fixed if there is no spin). 
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is a unitless expression of the difference between the bi- 
nary masses (which ranges between —1 and -1-1), and (j)r 
expresses the gravitational wave phase as a function of 
time. In addition to the variables S,i, and t/ir, the ex- 
pressions for the higher-order Hs also depend on rj. Note 
that the variables rj and S are not independent of each 
other: 



v = H^-s') 



(3.5) 



Since 77 does not uniquely specify the mass ratio (for ex- 
ample, 777,1 — b7n2 and 7712 = ^ttii yield the same 77), it is 



better to treat S a s independent and calculate rj from it 
using equation |3.5| . 

In the frame of the source, the emitted gravitational 
wave phase 4>s{ts) is related to the phase <j){t) of the bi- 
nary orbit as follows: 
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where loq is the orbital frequency at time t = and ts 
is coordinate time measured in the source frame. As 
discussed on page 579 in Blanchet, et. al. Q, the loga- 
rithmic correction term, while it is formally of 3/2 PN 
order, is actually of order 8/2 PN relative to the lead- 
ing term 4>(ts), and thus will be very small when e ^ 1. 
Therefore, we can safely make the approximation 



Mt) ~ (t>{t) 



(3.6b) 



which significantly simplifies the calculations. 

In the frame of the detector, the received wave is 
Doppler-shifted due to both the detector's orbital mo- 
tion and the source's radial motion with respect to the 
solar system, which we assumed was entirely cosmologi- 
cal. The received wave frequency ajj. = d(f>r/dt depends 
on the binary's orbital frequency ujs = d(j)/dts as follows: 
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Integrating both sides of this expression yields: 
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where 



/o(i) = / LUsSinint- <^)dt 



(3.8a) 
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and where z is the cosmological redshift factor, t is co- 
ordinate time in the detector frame (note that ts — 
[1 + z]~^t), 00 is the phase of the received wave at time 
t = [1 + z]ts — 0, (j>os is the phase of the binary's orbit at 
that instant, uJs = ^s{ts) is the angular frequency of the 
binary's orbit in its own frame at time ts, & and <& are 
the source's angular coordinates relative to the ecliptic, 
n is the angular frequency of the earth's orbit around the 
sun, and R is the radius of the earth's orbit. (We treat 
both H. and R as being constant; th us f lR is the Earth's 
constant orbital speed.) Equations 3.8 implicitly define 
<& = to be the same direction as that of a displace- 
ment from the sun of the detector array's center of mass 
at time t = 0. (As mentioned before, the polar angle Q 
corresponds to the zenith of the ecliptic.) Note also that 
since the earth's orbital speed fiR/c ~ 10~^ ^ 1, the 
Doppler shift has only been computed to first order in 
nR/c. 



According to Blanchet, et. al. [Q, the orbital phase, in 
turn, is given by: 



where 
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where tc is the time to coalescence from ts = 0, as mea- 
sured in the source frame. (Specifying tc is equivalent to 
specifying the initial separation of the orbiting masses.) 
Taking the time-derivative of this function gives the or- 
bital angular frequency cus = d(j)/dts in the source frame. 
Note that Us is specified in terms of time ts in the source 



frame, so in the integration described in equation 3.81 
the value of the integration variable t has to be trans- 
formed to ts in order to evaluate Wg correctly. 



B. Handling Cosmological Distances 

When dealing with cosmologically distant sources, it 
is more practical to specify the distance to the source 
in terms of the luminosity distance Rl = [L/iTrFob]^^'^ 
(where L is the source's intrinsic luminosity and Fob is 
the observed flux) or the redshift factor z instead of the 
geometrical distance r. In a flat, Friedmann (A = 0) 
universe, r, Rl, and z are related as follows: 

Rl 



1 



RlH = 2[l + z-VTT^ (3.10) 



where H is the current value of the Hubble constant (we 
used H — [14 X 10^ y]~^ in our calculation). We will con- 
sider the luminosity distance Rl to be the fundamental 
distance variable: if in some cases it is more practical to 
work with z, one can easily use equation 3.10 to change 
variables. 

Equations 3.1 to 3.10| therefore allow us to calculate 
(to 4/2 PN order) both the plus and cross polarizations 
of the received gravitational wave (even if the source is 
at cosmological distances) if we know the eight quantities 
T, tc, Rl, S, i, (po, 6, and <i>. 



C. The Detected Signal 

We work in the low-frequency limit in which the phase 
change S(j){t) of the detected laser tracking signal is di- 
rectly proportional to the amplitude of the wave: 
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h{t) ^ F+it)h+it) + F^{t)h^it) (3.11) 



where v is the laser frequency and T is the light-travel- 
time along the interferometer arm. The time-dependent 
functions F^{t) and Fx{t) are the beam-pattern func- 
tions for the interferometer pair. For both the precessing- 
plane case and the ecliptic plane case, the arms in a given 
interferometer pair make a 60° angle with respect to each 
other, and the beam-pattern functions for such a pair are 
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where Qd and (/)£> are the instantaneous angular coordi- 
nates of the source measured relative to the frame of the 
detector and V'u specifies the orientation of the binary's 
principal polarization axes around the fixed line of sight. 
The latter variable can be understood as follows. If the 
orbital inclination i to the line of sight is not zero or tt, the 
quasicircular binary orbit will look elliptical to a viewer 
at the detector. The angle f/'D specifies the orientation 
of the major axis of the ellipse as viewed by this observer 
measured in the plane perpendicular to the line of sight 
and from a reference direction that is perpendicular to 
both the line of sight and the normal to the plane of the 
detector array. If L, n, and zy) are unit vectors parallel to 
the binary's conserved angular momentum, the direction 
of the line of sight, and the normal to the detector array, 
respectively, then 
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In the ecliptic-plane case, zd is constant, so we will 
have ij^D = "0, which defines the fixed angle of the binary 
orbit's major axis relative to the plane of the ecliptic. In 
the precessing-plane case, z^ varies with time, so we will 
end up having to specify ■0D as a function of t, 8, $, and 
the fixed angle -0. 

In both the precessing-plane and ecliptic plane case, 
the equilateral arrangement of detector satellites means 
that the detector array possesses two independent pairs 
of interferometer arms, so we will have two independent 
signals hk{t) (where fc = 1, 2). Because one pair of arms 
is rotated 60° relative to the other in the plane of the 
detector array, the value of (pD will depend on the choice 
of interferometer pair as well, and so is denoted 4>D,k- 



For the ecliptic plane case, the zenith of the detector 
plane is the same as the zenith of the ecliptic, so O^, = 
<d and, as we said, -0-0 = "0- If the satellites orbit the 
Earth with an angular frequency of lj^ , then the apparent 
azimuth of the source relative to the detector arm will be 
given by 

0_D,fc = a.k{t) where ak{t) = $ - ujdt + aofc (3.14) 

Here the constant aofe specifies the orientation of the 
interferometer pair at i = 0. While we must have 
ao2 = oiQi + -I, the constant agi = Q^o is arbitrary. In 
terms of these variables, the beam-pattern functions be- 
come 

F+,fc = ^ [i(H-cos2e)cos2afcCOs2-0 

— cos 9 sin 2Q;fe sin 2-0] (3.15a) 

Fx ,fe = ^ [i (l-f cos^ 9) COS 2ak sin 21/- 

+ cos9sin2afc cos20] (3.15b) 

In the precessing-plane case, the orientation of the de- 
tector plane changes with time, making the expressions 
for the beam pattern functions more complicated. As 
shown in Cutler ^ (and verified by us), to get the beam- 
pattern functions for the prec essing plane case, we must 
substitute into equations 3.12 the quantities 
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where in this case 

akit) = nt + aQk (3.18) 

and /3{t) specifies the angular position of the detector 
array's center of mass in the plane of the ecliptic. Tak- 
ing account of the way that we have defined <& = 0, the 
quantity f3{t) is given by: 

P(t)=nt~<P (3.19) 

One can also show that i/i/^ in this case is given by 
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where 
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To summarize, we see that the detected signal from a 
inspiralling binary source depends on nine parameters r, 
tc, Rl (or z), rj (or S), i, 0Oj V": ©i and $, each of which 
we would like to determine from the observed signal. In 
addition, we have one parameter ag that we will know in- 
dependently of the signal (once the array is launched, we 
should know the orientation of the interferometer arms 
as a function of time and thus ao at the given instant we 
define to be < = 0). 



D. Statistical Methods 

The goal of our analysis is to estimate the accuracy 
with which the parameters that determine h{t), espe- 
cially the positional parameters & and <I>, may be de- 
termined from the data. The method we will use is 
that of linear least-squares parameter estimation. In this 
method, we assume that approximate initial values have 
been found for the m parameters q^ that determine the 
signals. We also assume that n observations hi of the 
waveform have been made and that the difference be- 
tween the observed waveform and the waveform predicted 
by the current approximate values of the parameters is 
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where v^ is the noise associated with each measurement 
and Ago is the error in the current value of the parameter 
qa ■ The goal of least squares parameter estimation is to 
find the adjustments Aqa to the parameters that will 
minimize the Vi errors in a least-squares sense. If we 
define 



Xa = Aqa and G^a = T—^ 
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then the total squared error (also called the variance) 
may be written 



V = ViUi = (yi - GiaXa){Vi - GicXc) 
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where sums over i from 1 to n and over a and c from 
1 to m are implied by the repetition of the subscripts. 
The least squares requirement of minimizing the variance 
proceeds by requiring that 

^ = -2y^Ga + 2G^bG^aXa = (3.25) 

dxb 

If we define an to x tti square information matrix 

■A-ab — GiaGib (3.26) 

then the information equation \i.2^ may be written 

Aab'Xa = ViGib (3.27) 

giving a solution for Xa of 

Xa = Ajy,G,b (3.28) 

where A~i^ is the matrix inverse of Aab ■ 

The uncertainty in Xa that results from such a proce- 
dure is related to the uncertainty ay in the observations 
yi. It is found by defining the covariance matrix as the 
expectation value of the mean-squared deviation of Xa 
from its expected value: 



Cab = {[Xa - {Xa)][xb - (Xb)]) 



(3.29) 



where (...) represents the expectation value of the quan- 
tity inside the brackets. Since only Xa and yi are random 
variables, the expectation value of Xa is given by 

(xa) = A-^G,, (y,) (3.30) 

Combining equations 3.28|, 3.30| and ^.2g|, we get 



Gab ^ A-^G.,A-^G,d {[y^ - {y^)] [y, - (%)]) (3.31) 

If the errors in the observations arc un iformly distributed 
and uncorrelated (see subsection [II F ) , then the expecta- 
tion value of the mean squared error in the observations 

is 



{[vi~ {yi)][y3- {V])]) = Sij(jl 



(3.32) 



where Sij is the Kronecker delta. Substituting equation 
3.32 into equation 3.31 , and remembering that AcdA'i^^ = 
6cbi we can write 



Gab — l^y^ab 



(3.33) 



According to equation 3.29 the expected errors a a in the 



parameter values are given by the the diagonal elements 
of the covariance matrix Gab'- 



\J{[Xa- {Xa)][Xa " {Xa)]) 



(3.34) 



Therefore, equation p. 33 implies that the expected errors 
are equal to 



(3.35) 



Ca ~ CTyV Aaa 

(no sum over the index a is implied here). 

E. Calculating the Derivatives 



In order to construct the information matrix for the 
uncertainty calculation, we need to calculate the par- 
tial derivatives of the hypothetical observed signal hk (t) 
with respect to each of the nine unknown parameters 
T,tc,RL,S,i,(j)o,''P,Q, and $. We evaluated all of these 
derivatives analytically to get specific expressions for 
dhi;/dq (where q is any one of the nine variables) as 
a function of time. This is mostly a matter of sim- 
ple (though tedious) differential calculus. The following 
notes give an overview of the calculations. 

Consider the derivative with respect to an arbitrary 



variable q. According to equation 3.11 



dhk 
dq 



dK 



+,fe 
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dq 
dq 



■F. 



hx+F: 



dh+ 
dq 



x,fe' 



dq 



(3.36) 



Of the partials of the beam-pattern functions F, only 
dF/dQ, dF/d^, and dF/dip are nonzero, so evaluating 



these partials is relatively straightforward. If we look at 
equation 3.1, we see that 



dhn 



dq dq 



-2 u(0) 



-3fr(l/2) 



6 Er(4/2) 
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+ ,x 

de 
dq 

(3.37) 



where ho = 2Tcr]/5r. 

Now, many of the quantities appearing in this equation 
are written in terms of the variable rj, but the more fun- 
damental variable describing the difference between the 
masses involved in the binary is S. Since ry is a simple 
function of the fundamental variable S alone, we can use 
the chain rule to evaluate partial derivatives with respect 
to S as follows: 



a/(77, ...) dfiv,...)d7j df{r,,...) 
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drj dS 



dr] 



(3.38) 



Similarly, note that ho is stated in terms of the geometric 
distance to the source r, but we are taking the fundamen- 
tal distance variable to be the luminosity distance Rl , so 
we evaluate the partial derivative of ho with respect to 
Rl as follows: 



dho 
dRr 



dho dr dho 1 



dr dRL 9r 1 -I- z 



(3.39) 



The expansion parameter e explicitly depends on t, 
and (through the angular frequency lo) on the variables 
tc and rj (or 6), so the only partials we need to evaluate 
the second line of equation 3.37 are: 
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(3.40a) 



(3.40b) 



(3.40c) 



the partials of the Hs. These functions, which in the 
higher-order cases can become quite complex, depend 
only on i and d and/or ry explicitly, but also depend on the 
received phase (j)r, which in turn depends on r, ic, (5, (/loj © 
and <&. Therefore, generally we have to evaluate, for each 
of the ten H terms. 



dH dH ^ 9H _ dH d<j)r 
di ' dS ' dq d(j)r dq 



(3.41) 



Evaluating the derivatives of (jjr is complicated by the 
fact that one cannot easily evaluate the Doppler-shift in- 
tegral in equations 3^ analytically because of the com- 
plicated time-dependence of the source-frame orbital fre- 
quency LOs- Therefore, we expressed each of the partial 
derivatives of (pr in terms of an integral equation: for 
example 



d(j)r 

~d7 



1 



90 

a7 
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sin 8 



~dT 



su\{nt - $) dt 



(3.42a) 
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l + z 



nR 

c 



cos Q ujg sm{flt — <i>) dt 
Jo 



(3.42b) 



We then evaluated these integrals at each time step using 
a fourth-order Runge-Kutta integration scheme. Since 
ujs is a fairly slowly-varying function of time (except in 
the final stages of coalescence, where the whole post- 
Newtonian power-series approximation breaks down any- 
way), the Runge-Kutta approach should yield excellent 
estimates for the values of these integrals. 



F. Noise Curves 



Ultimately, the detector's position sensitivity is deter- 
mined by the ay noise in the detector. We took the total 
noise for the covariance study to be 



al - SM) A/ 



(3.43) 



where the bandwidth A/ — l/2dt {dt being the time be- 
tween samples), and 5'„(/) is the spectral noise density. 
The total rms noise in the detectors, sampled at intervals 
dt, would of course be the integral of 5'„(/) over the band- 
width. By taking instead the simple product of A/ and 
Sn (/) at the primary signal frequency, we account for the 
temporal correlations created in the sampled data by the 
larger low-frequency noise components (see Appendix for 
details). 

To estimate the spectral noise density >5'„(/), we have 
taken the published noise curves for the proposed space 
missions. For the ecliptic-plane scenario, we used the 
noise curve for OMEGA; for the precessing-plane case, 
we use the noise curve for LISA. The instrumental noise 
curves for both missons are dominated by two sorts of 
noise, acceleration noise and position noise. The genera- 
tion of sensitivity curves with these types of noise present 
is discussed in Larson, Hiscock, and Hellings p]. In addi- 
tion, all instruments in the low-frequency band will have 
to contend with a background of clutter formed by un- 
resolvable close compact binaries in the Galaxy g. The 
total spectral noise density function is thus given by 



for q = T, tc, S, (j)o, 6, and $. 



Sn{f) = [SaU)r^ + 5,(/)](l + 27r/r)2 + Sc{f) (3.44) 



where Sa{f) is the strain noise produced by accelera- 
tion noise, Sx{f) is the strain noise produced by position 
noise, and T is the hght-travel time between detectors. 
The (1 + 27r/T)^ factor arises because the sensitivity of 
the instrument to gravitational waves is reduced when 
the gravitational wave wavelength is less than the base- 
line distanc e bet ween spacecraft. The values of the terms 
in equation p. 44 are given in Table |[ In addition, in both 
missions we are assuming that we are looking at clutter 
noise of 



Scif) 



(2.07 X 10- 
(4.73 X 10 
(1.41 X 10 



^43 HzO-9)/- 

-61 Hz6-5)/- 
-47 Hzl-6^ 



')r 



1.9 
7.5 
2.6 



/</l 

h </</2 

/>/2 



(3.45) 

where /i = 7.1 x 10"^ Hz and /a = 1.8 x 10"^ Hz. The 
total spectral noise curves resulting from equations 3.44 
and 3.45 are shown in Figure 0. 



G. Summary of the Algorithm 

To perform the actual calculations, we wrote a FOR- 
TRAN program that does the following: 

1. The program reads values for the fundamental pa- 
rameters TOi + m2,tc,RL (or z), 6,i,(f)o,ip,Q, and 
$ for a hypothetical source from an input file. In 
addition, the input file also specifies the size of the 
time step dt to be used and the detector array's 
initial orientation angle ag. (We generally used dt 
= 60 s in our runs so that dt was somewhat smaller 
than the shortest period for waves that either de- 
tector can easily register.) 

2. The program then uses the specified values to set up 
a number of constants used in the calculation, and 
initializes the Runge-Kutta scheme used to evaluate 
the integrals in the expressions for (f>r and its partial 
derivatives. 

3. It then executes a loop that at each time step (in- 
cluding i = 0) 

• calculates the noise level 5'„ at the current ob- 
served gravity wave frequency, 

• calculates values of 0r and its partial deriva- 
tives, 

• calculates hk{t) and its partial derivatives for 
each detector, 



• adds the information from the current time 



step to matrix A*^ — Ag b/al (where Aab is as 



defined in equation 3.26| ) as follows: 

A:,{t) = A:,{t-dt) 



dt 



Snif) , 



V ^^ (3.46) 
^ , dqa dqb 



where qi corresponds to the z-th fundamental 
parameter. 

• evaluates the Doppler-shift integrals needed 
for the next time step. 

The loop ends when the time to coalescence is 
smaller than 10 times the current orbital period or 
after one full year of observation, whichever occurs 
sooner. 

4. The program then uses a standard matrix-inversion 
package to invert the matrix y4*j, to get Cab = 
'^y^ab ~ ( ^*h)~ ^ (see equation 3.33). According 
to equation |3.35 , the diagonal elements of this in- 
verted matrix correspond to the uncertainties in the 
values of the fundamental parameters that we could 
determine from the wave train generated by the hy- 
pothetical source in question. If the matrix cannot 
be inverted, the program eliminates the parame- 
ter causing the problem and computes uncertainties 
for the remaining parameters. Where possible, the 
program simultaneously estimates uncertainties for 
all nine of the fundamental parameters, not just 
the angular position parameters. 

5. It then displays the calculated results for the un- 
certainties, including the solid-angle uncertainty in 
the angular position of the source. 

Running this program repeatedly with different param- 
eter values allows us to learn about how the parameter 
uncertainties might depend on the characteristics of the 
source. 

We tested the program by making certain that in the 
low-mass limit, its results agreed completely with a much 
simpler program we had written earlier to handle the case 
of strictly monochromatic waves from low-mass binaries. 
Even though significant differences between our assump- 
tions and those of Cutler R] and Cutler and Vecchio pP] 
made exact comparisons difficult, we also tested our pro- 
gram by comparing its results to their published results 
for the precessing-plane case, and found agreement typi- 
cally within a factor of three. 



IV. RESULTS 

The results in principle depend on ten parameters: the 
nine fundamental (unknown) parameters and the pre- 
sumably known initial orientation parameter ao- Ob- 
viously, it is difficult to sample such a large parameter 
space and display its results meaningfully. We did find 
that in spot checks the uncertainty in the source's an- 
gular position seemed to be generally insensitive to the 
parameters (f>o, "0? Q;oj and $, and that the results were 
basically the same for polar angles O above and below 
the ecliptic that made the same angle with the ecliptic. 
In this paper, we will focus primarily how the angular 



position uncertainty depends on the total mass of the 
source, or its frequency, and the polar angle Q. 

Figure ^ shows the base-ten logarithm of the angular 
position uncertainty (expressed in stcradians) as a func- 
tion of for essentially monochromatic compact binaries 
at selected frequencies ranging from 0.1 mHz to 100 mHz. 
The masses were chosen so that tc was long compared to 
a year, the distance was chosen so that the signal-to- 
noise ratio in the ecliptic case was roughly 10, and the 
inclination was chosen so that cosi was about 0.8. 

Note that in the ecliptic-plane case (solid curves), there 
is a monotonic improvement in angular position resolu- 
tion with increasing frequency, as the Doppler shift be- 
comes easier to track. Indeed, in the strictly monochro- 
matic limit, the angular position uncertainty should fall 
strictly as f~^ (see page 99 of H), and our results are 
consistent with this. Note also that the angular position 
uncertainty increases as one approaches 8 = 90° (that is, 
the plane of the ecliptic) . This is because in the ecliptic 
plane case, the only source of in forma tion for Q comes via 
the Doppler shift (see equation 3.8a). Since this effect is 



proportional to sin Q the partial derivative with respect 
to 8 is proportional to cos 8 which goes to zero at 8 = 
90°. There is thus no 8 information in the information 
matrix at 8 = 90°. 

In the precessing-plane case (the dotted curves) , there 
is moderately good angular resolution even at low fre- 
quency because of the extra information provided by the 
precession of the array plane. The angular resolution im- 
proves as frequency increases because of the additional 
information provided by the Doppler shift. The angular 
resolution in this case proves to be relatively independent 
of 8 (the precession of the array plane means that the 
ecliptic plane is not as special as it is in the ecliptic case). 

Figure shows the angular uncertainty as a function of 
8 for various large-mass coalescing sources. All of these 
sources have equal-mass partners (at 10^, 10^, 10^, and 
lO"' solar masses respectively), are located at redshift z 
= 1, have an inclination whose cosine is 0.8, and are ob- 
served for the last year before coalescence. Note that in 
the precessing-plane case (the dotted curves) the angu- 
lar resolution is fairly good even for the lowest-frequency 
cases (the cases with the largest masses) . The resolution 
stays relatively constant as the mass decreases, because 
even as the signal strength drops, the noise curve also 
drops. Again the angular position uncertainty is rela- 
tively independent of 8. 

In the ecliptic-plane case, the angular resolution again 
improves as the frequency increases (as expected), but 
curiously, and in sharp contrast to the monochromatic 
case, the angular resolution improves significantly as 8 
approaches 90°. We will discuss this further in the next 
section. 

Figure [^ shows graphs for middle-mass black holes 
falling into supermassive black holes. One can see that 
the results are qualitatively similar to the equal-mass 
black-hole mergers shown in Figure m 



V. DISCUSSION OF THE RESULTS 

In the ecliptic plane case, the angular position uncer- 
tainty decreases significantly as we approach the eclip- 
tic plane in the massive black hole cases shown in Fig- 
ure 0, but the uncertainty increases significantly in the 
monochromatic case shown in Figure |3| What could be 
so different about the two cases? It turns out that the 
crucial difference in determining this effect is that in the 
massive black hole case, the higher-order post-Newtonian 
corrections to the waveform amplitude are significant 
(because the expansion parameter s in the power series 
is large), while in the monochromatic case, they are not. 
Figure clearly shows that when we artificially turn off 
the higher-order terms, the angular position uncertainty 
increases with 8 for the massive black-hole case, just 
as it does for the monochromatic case. Our runs show 
that virtually all of the difference comes from the lowest- 
order nonzero harmonic above the fundamental: adding 
only this term to the fundamental produces essentially 
the same curve we get when we include all higher-order 
terms. 

Why is this? We may consider the waveforms seen at 
each of the two detectors as depending on three effects. 
The first effect is the monotonic incr ease in the frequency 
of the source, as given by equations 3^. If this increase 
is written by expanding the frequency in a Taylor se- 
ries, then the behavior can be expressed in terms of the 
derivatives wq, too, '^o, etc. As seen in equations 3.9, these 



derivatives are linked to the basic variables rj, r, and tc- 
Observation of the time series will thus determine the fre- 
quency derivatives and thus rj, r, and tc independently of 
any other features of the waveforms. The second effect 
is the variation of the waveform with orientation of the 



detector, as given by equations 3.15 . These form factors 
depend explicitly on 8, $ (through a^), and "0. The 
form factors may not be seen directly in the waveform, 
but only in convolution with the third effect, the ampli- 
tudes of the two polarizations of the waves, as given by 
These /i-functions are determined by the 



3.1 



equation 

already- known r, 77, and w, and by the unknown parame- 
ters r, i, and (j)o. There are thus six unknown parameters 
that must be determined from the time series alone, with- 
out any help from the change in frequency with time. If 
only the fundamental fr equen cy o f the gravitational wave 
were present (equations p. 3a and 3.31 ) , then each detec- 
tor would see only a single gravitational wave frequency 
whose amplitude and phase would be the only observ- 
ables. For two detectors, there would be two amplitude 
observables and two phase observables, but this would 
not be enough to determine the six unknown parameters, 
8, $, ijj, r, i, and (pQ. However, if a second harmonic of 
the wave is in cluded (the harmonics ij'^/-^) of equations 
3.3c and 3.3d ), then a Fourier analysis of the detected 
signals would determine amplitudes and phases for both 
harmonics. Since the mix of phases and amplitudes be- 
tween the two harmonics depends on i and (jjo, there will 
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be nontrivial information in these additional amplitudes 
and phases, and all six gravitational wave parameters can 
be determined from these eight signal response parame- 
ters. The ability of the higher harmonics to add informa- 
tion, however, is quenched as -^^ 0, because the form 
factors depend on 8 only through cos O and also because 
$ and if] become degenerate near the pole of the ecliptic. 

In any case, Figure |6| makes it clear that including the 
higher-order terms makes a huge difference for values of 
Q close to the ecliptic plane. It is worth noting in this 
context that if massive-black hole sources are distributed 
evenly about the sky, half will lie within ±30° of the 
ecliptic plane, where the effect of including the higher- 
order terms is most significant. 

This result has not been previously reported in the 
literature. Cutler m, and Cutler and Vecchio p] only 
model the gravitational waveform itself to lowest post- 
Newtonian order, leading them to conclude (appropri- 
ately, considering the limitations of that model) that 
a detector array in ecliptic-plane configuration "has es- 
sentially no angular and distance resolution" for mas- 
sive black hole mergers (see page 105 of 0). However, 
the results reported here indicate that we need not be 
so pessimistic: indeed, for Q very near the ecliptic, the 
angular resolutions in the ecliptic-plane and precessing- 
plane cases are quite comparable as long as we include 
the higher-order terms. 

It is interesting to note that we do not observe any- 
thing comparable for the precessing-plane case: we find 
that artificially suppressing the higher-order terms in the 
waveform does not change the angular uncertainties very 
much. In the precessing plane case, the fact that the nor- 
mal to the detector plane moves around the sky means 
that 9d and ipo depend on time in a characteristic way 
even though Q and ij) are fixed. This provides informa- 
tion that allows one to disentangle the effects of the latter 
parameters from those of i and Rl even without the ex- 
tra information provided by the higher-order terms in the 
wavefunction. 



VI. CONCLUSIONS 

While in the main the results from this study validate 
much the work that has already been done on the an- 
gular position resolution of space-based detectors, this 
work underlines the unexpected but critically important 
role that higher-order post-Newtonian corrections play 
in determining parameters in massive black-hole mergers, 
particularly for a detector array in the ecliptic-plane con- 
figuration. This suggests that future research into how 
one determines parameters of massive black-hole mergers 
from an observed wave train should definitely take ac- 
count of these higher-order post-Newtonian corrections, 
at least if one is considering an ecliptic-plane detector 
configuration. 

It is also worth noting that our results indicate that 



neither detector configuration is able to yield an angular 
position uncertainty much less than about 10~^ stera- 
dians for massive black-hole mergers (see Figures and 
H). For comparison, the full moon occupies a solid an- 
gle of about 6 X 10"^ sr. This means that even though 
either kind of detector can easily register massive-black- 
hole mergers at cosmological distances, neither will be 
able to locate the source precisely enough (by several 
orders of magnitude) to indicate the galaxy where the 
merger is taking place. One can therefore only locate the 
merger if one can also observe the electromagnetic sig- 
nature of such an event. To our knowledge, little is yet 
definitively known about what the electromagnetic sig- 
nature of a merger might look like. Using the position 
information provided by either kind of detector in prac- 
tice (or even knowing whether the information has any 
practical use at all) will thus require more research into 
the astrophysics of black-hole mergers. 
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APPENDIX: HANDLING RED NOISE 

The gravitational wave sources we have considered in 
this paper are all found at frequencies on the left of Fig- 
ure 2, where the noise spectrum is falling roughly like 
f^^. Noise such as this, where there is more power 
at low Fourier frequencies than at high frequencies, is 
termed "red noise". It is generally known that param- 
eter estimation, including spectral estimation, produces 
biased estimates when the background noise is not white. 
[0,|| The solution to this bias problem lies in a pro- 
cess called "prewhitening" . In this process, the data are 
passed through a filter that generates a new time series 
whose noise spectrum is flat, and the parameter estima- 
tion is performed using this filtered data. The estimates 
will then be bias-free. In this appendix, we show how 
to calculate the signal-to-noise ratio (SNR) for such a 
process, assuming the signal is dominated by a single 
frequency. 

Let us consider a time series y{t) — h{t) + n{t), where 
h{t) is the signal and n{t) is the red noise. To prewhiten 
the data, the time series is passed through a liner filter 
to produce: 



x{t) = F{t) * yit) 



loFit- 



T)y{T)dT 



(Al) 



We write the effect of this filter on the signal h{t) as 
g{t) — F{t) * h{t) and on the noise as m{t) = F{t) * n{t). 
The filter F is chosen so that the power spectrum of ■m{t) 
will be fiat. Since the filter is a convolution integral, the 
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effect of the filter in frequency space yields simple prod- 
nets, g{f) = F{f)h{f) and to(/) = F{f)n{f), where 
F{f) is called the transfer function of the filter. [For 
example, if the power spectrum of the detector noise 
Snif) = "-(Z)^ were exactly proportional to /~^, then 
the linear filter would be F(t) ~ d?/dP whose transfer 
function is F{f) — 47r^/^. The power spectrum of m{t) 
would therefore be fiat.] 

The SNR for the filtered data is given by (SNR)^ = 
{g"^ {t)) / {w? {t)) , where the angular brackets denote a 
time average. The mean squared signal and noise 
strengths may be written in terms of their spectral den- 
sities to give 



(SNR) 2 



I}"S,if)df 

Jj^ F^if)S,if)df 

SmifH — Jl) 



(A2) 



where, in the last step, the fact that Sm{f) — const has 
been used to complete the integral in the denominator 
and the fact that g{f) — F{f) * h{f) has been used to 
expand the numerator. 

Now let us assume that the signal h{t) is essentially 
monochromatic, at frequency /o, so that in the numera- 
tor we will have Sh{f) = S{f—fo)Sh{fo)- We also assume 
that the data have been previously high-pass filtered with 
cutoff at Jl = fo- In the denominator, because Sm{f) 
is constant, its relation to S'„(/) may be worked out at 
any frequency desired. If we choose /o as that frequency, 
then we have 



(SNR)2 



FHfo) li" S{f - fo)S,ifo)df 



FHMSnifoKfH 



fo) 



5„(/o)(l/2At) 



(A3) 



where, in the last step, we have recognized {h{t)'^) as the 
integral of the spectral density Sh and fn ~ l/2Ai as 
the Nyquist frequency, assumed to be much higher than 
fo. This is the form ula w e use for noise in our detector, 
as given in equation 3.4E . 
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FIG. 1. Spacecraft configurations for the ecliptic-plane and 
precessing-plane options. Figure nla shows top and side views 
of the three detector spacecraft in the ecliptic-plane option. 
The arrow shows the normal vector to the plane of the orbit, 
which is fixed in this case. Figure wo shows top and side 
views of the three detector spacecraft in the precessing-plane 
option. In this case, the spacecraft do not orbit the earth but 
rather orbit the sun in such a way that the satellites form 
an apparently rigid equilateral triangle that rotates clockwise 
around the array's center, which in turn orbits the sun. The 
plane of the array is tilted 60° from the plane of the ecliptic 
(which is why the circular trajectories of the satellites look 
elliptical in the top view). The projection on the ecliptic of 
the array's normal vector always points toward the sun. 



FIG. 2. Graph of the noise spectral density S„{f) that was 
used in the covarance analysis. The solid curve is for the 
ecliptic-plane case (using OMEGA noise estimates) while the 
dashed curve is for the precessing-plane case (using LISA noise 
estimates) . 



FIG. 3. Angular uncertainty as a function of B for a com- 
pact binary emitting gravitational waves with various differ- 
ent (but essentially constant) frequencies. The stellar masses 
were chosen in each case so that tc was much longer than one 
year, and the distance was chosen so that the signal-to-noise 
ratio was about 10 in the ecliptic-plane case (the signal am- 
plitudes are specified in each picture). The other parameters 
were set somewhat arbitrarily as follows: qq = 0, i = 36.9° 
(cosi = 0.8), <!> = 268.5° = 4.69 rad, ip = 114.6° = 2.0 rad. 



FIG. 4. Angular uncertainty as a function of O for various 
equal-mass pairs of supermassive black holes undergoing co- 
alescence. In each case, the coalescing pair was assumed to 
be at redshift 2 = 1 and the time to coalescence tc was cho- 
sen to be one year. The initial gravitational wave frequency 
shown in each graph is the frequency at the beginning of the 
one-year observational period, and the final frequency is the 
frequency when the calculation ended roughly 10 orbits be- 
fore coalescence. The other parameters were set as in Figure 
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FIG. 5. Angular uncertainty as a function of O for two 
situations where a medium-mass black hole coalesces with a 
supermassive black hole. Other parameters were set as spec- 
ified in Figure H. 



FIG. 6. Angular uncertainty as a function of for coa- 
lescing massive black holes, each with a mass of 10*' solar 
masses (this is the same situation shown in Figure w>)- This 
illustrates that the higher-order terms in the post-newtonian 
approximation for the gravitational waveform have a large ef- 
fect on the angular resolution in the ecliptic-plane for values of 
Q near the ecliptic plane: artificially suppressing these terms 
leads to very poor resolution near the ecliptic plane, while in- 
cluding them leads to an angular resolution that is nearly the 
same as for the precessing-plane case when O is close to the 
ecliptic. It turns out that virtually all of the difference is due 
to the lowest-order nonzero harmonic term above the funda- 
mental: adding this term alone produces a curve essentially 
identical to the lower curve shown here. 



TABLE I. Table showin g the magnitudes of the spectral 
noise constants in equation 3.44 for the OMEGA and current 
LISA detector configurations. 



Term 



OMEGA 



LISA 



'SJJ) 5.8 X 10""' Hz'^ 

&(/) 1.26 X lO"*! Hz-i 

2tvT 21 s 



2.3 X 10""^ Hz'' 
1.26 X lO-'i Hz" 
100 s 
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